
# ____________________________________________________________
# Replication for Paper Results
# ____________________________________________________________

rm(list = ls())
source("code/_preamble.R"); preamble()
source("code/_functions.R")

kompally_dta <- read_csv("data/election/kompally.csv")

# Table 1 --------------------------------------------------------------------
# (Balance Table)
sims <- 100000
sims <- 10
# List variables to check balance on
covars <- c("total_eligible","male_eligible","female_eligible", "n_party", 
            "percentmuslim_method1", "percentmuslim_method2", "women_res", 
            "stsc_res", "bc_res", "turnout_ge")

# Create a balance table shell
balance <- (matrix(nrow=length(covars), ncol=4))
balance <- as.data.frame(balance)
colnames(balance) <- c("No FRT Means","FRT Means", "Difference in Means", "Exact P-Value")
rownames(balance) <- covars

# Loop to create balance table
set.seed(1234)
for(i in 1:length(covars)){
  ca <- kompally_dta %>% filter(facial_recognition==0) %>% subset(select = covars[i]) %>% as.matrix() %>% c()
  ta <- kompally_dta %>% filter(facial_recognition==1) %>% subset(select = covars[i]) %>% as.matrix() %>% c()
  test_ct <- t.test(ca, ta, na.rm = TRUE)
  balance[covars[i],"No FRT Means"] <- test_ct$estimate[1]%>%round(.,digit = 2)
  balance[covars[i],"FRT Means"] <- test_ct$estimate[2]%>%round(.,digit = 2)
  balance[covars[i],"Difference in Means"] <-  test_ct$estimate[1]%>%round(.,digit = 2) - test_ct$estimate[2]%>%round(.,digit = 2)
  dv <- kompally_dta %>% select(covars[i]) %>% dplyr::rename(dv = covars[i])
  kompally_dta <- cbind(kompally_dta, dv)
  
  # use randomization inference for the p-values for differences
  fisher <- conduct_ri(
    formula = dv ~ facial_recognition,
    declaration = declare_ra(N = 36, m = 10),
    assignment = "facial_recognition",
    sharp_hypothesis = 0,
    data = kompally_dta,
    sims = sims
  )
  balance[covars[i],"Exact P-Value"] <- summary(fisher)[3]
  kompally_dta <- kompally_dta %>% select(-c(dv))
}

# an error will be produced beause there is no exact p-value for turnout_ge 
# this is calculated separately below

# calculate randomization inference for turnout_ge (since it has less observations)
turnout_ge_ri <- conduct_ri(
  formula = turnout_ge ~ facial_recognition,
  declaration = declare_ra(N = 26, m = 7),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta %>% filter(!is.na(turnout_ge)),
  sims = sims
)

beepr::beep(2)

# put exact p-value into balance table
balance[10,4] <- summary(turnout_ge_ri)[3]

# Saving table
# balance %>% 
#   as_tibble() %>% 
#   mutate(`Exact P-Value` = round(`Exact P-Value`, 2)) %>% 
#   mutate(Variable = c("Total Eligible Voters", "Male Eligible Voters", "Female Eligible Voters",
#                       "Number of Parties", "Percent Muslim (1)", "Percent Muslim (2)",
#                       "Gender Reservation", "SC / ST Reservation", "BC Reservation",
#                       "General Election Turnout")) %>%
#   select(Variable, everything()) %>% 
#   xtable(digits = 2) %>% 
#   Kable() %>% KableSave(fname = "balancetable")

# Table 2 --------------------------------------------------------------------
# (Main Effect Table)

# Models
t1m1 <- felm(turnout ~ facial_recognition | 0 | 0 | 0, 
             data = kompally_dta)
t1m1 %>% summary(., robust = T)

t1m2 <- felm(turnout ~ facial_recognition + n_party + any_res | 0 | 0 | 0, 
             data = kompally_dta)
t1m2 %>% summary(., robust = T)

t1m3 <- felm(turnout ~ facial_recognition + pfemale + total_eligible + percentmuslim_method2 | 0 | 0 | 0, 
             data = kompally_dta)
t1m3 %>% summary(., robust = T)

t1m4 <- felm(turnout ~ facial_recognition + pfemale + total_eligible + 
               n_party + any_res + percentmuslim_method2 | 0 | 0 | 0, 
             data = kompally_dta)
t1m4 %>% summary(., robust = T)

# Randomization Inference for Table 2 Models 1-4 to get exact p-values
set.seed(1234)

t1m1_p <- conduct_ri(
  formula = turnout ~ facial_recognition,
  declaration = declare_ra(N = 36, m = 10),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)

t1m2_p <- conduct_ri(
  formula = turnout ~ facial_recognition + n_party + any_res,
  declaration = declare_ra(N = 36, m = 10),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)

t1m3_p <- conduct_ri(
  formula = turnout ~ facial_recognition + pfemale + total_eligible + percentmuslim_method2,
  declaration = declare_ra(N = 36, m = 10),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)

t1m4_p <- conduct_ri(
  formula = turnout ~ facial_recognition + pfemale + total_eligible + 
    n_party + any_res + percentmuslim_method2,
  declaration = declare_ra(N = 36, m = 10),
  assignment = "facial_recognition",
  sharp_hypothesis = 0,
  data = kompally_dta,
  sims = sims
)

# Save p-values
t1m1_pval <- summary(t1m1_p)$two_tailed_p_value %>% round(2)
t1m2_pval <- summary(t1m2_p)$two_tailed_p_value %>% round(2)
t1m3_pval <- summary(t1m3_p)$two_tailed_p_value %>% round(2)
t1m4_pval <- summary(t1m4_p)$two_tailed_p_value %>% round(2)

# Saving table
# stargazer(t1m1, t1m2, t1m3, t1m4,
#           keep = c("facial_recognition"),
#           covariate.labels = c("Facial Recognition Technology"),
#           se = list(t1m1$rse ,t1m2$rse, t1m3$rse, t1m4$rse),
#           star.cutoffs = NA,
#           digits = 2,
#           # ci = T,
#           # p = list(c(t1m1_pval, t1m2_pval, t1m3_pval, t1m4_pval)),
#           report=('vc*s*p'),
#           dep.var.labels.include = FALSE,
#           keep.stat = c("n"),
#           title = "Effect of Facial Recognition Software on Voter Turnout",
#           add.lines = list(c("Mean of DV", rep(round(mean(kompally_dta$turnout), 2), 4)),
#                            c("Election Controls", "N", "Y", "N", "Y"), 
#                            c("Polling Station Controls", "N", "N", "Y", "Y")),
#           float = F,
#           table.placement = "H",
#           omit.table.layout = "n",
#           notes.append = F,
#           no.space = T,
#           label = "table:maineffect",
#           out = "tables/maineffect.tex") 
# 
# # Replace regular p-values with exact p-values (Fisher)
# ri_pval = paste0("& [", t1m1_pval, "]", "$^{*}$", " & [", t1m2_pval, "]", "$^{***}$",
#                  " & [", t1m3_pval, "]", " & [", t1m4_pval, "]", "$^{**}$", "\\")
# 

# Figure 1 --------------------------------------------------------------------
# (Marginal Effects Plot)

d <- cbind.data.frame(Y = kompally_dta$turnout,
                      D = kompally_dta$facial_recognition,
                      X = kompally_dta$percentmuslim_method2,
                      Z1 = kompally_dta$pfemale, 
                      Z2 = kompally_dta$total_eligible, 
                      Z3 = kompally_dta$n_party, 
                      Z4 = kompally_dta$any_res
)

m_no_control <- interflex::interflex(Y = "Y", D = "D", X = "X", 
                                     data = d, estimator = 'binning', nbins = 3, bin.labs = FALSE,
                                     vcov.type = "robust", 
                                     theme.bw = T,
                                     Ylabel = "Turnout",
                                     Xlabel = "Percent Muslim",
                                     Dlabel = "FRT",
                                     cex.lab = 1,
                                     cex.main = 1,
                                     ylim = c(-65, 10),
                                     main = "Model without control variables")

m_controls <- interflex::interflex(Y = "Y", D = "D", X = "X", 
                                   Z = c("Z1", "Z2", "Z3", "Z4"),
                                   data = d, estimator = 'binning', nbins = 3, bin.labs = FALSE,
                                   vcov.type = "robust", 
                                   theme.bw = T,
                                   Ylabel = "Turnout",
                                   Xlabel = "Percent Muslim",
                                   Dlabel = "FRR",
                                   cex.lab = 1,
                                   cex.main = 1,
                                   ylim = c(-65, 10),
                                   main = "Model with control variables") 
# Saving
m_no_control$figure + m_controls$figure
ggsave("figures/fg1.pdf",
       height = 5,
       width = 10)

